home *** CD-ROM | disk | FTP | other *** search
/ Enter 2006 September / Enter 09 2006.iso / Internet / SpamExperts Home 1.1 / SpamExperts Home.exe / lib / spamexperts.modules / spambayes / chi2.pyc (.txt) < prev    next >
Encoding:
Python Compiled Bytecode  |  2006-07-14  |  5.6 KB  |  202 lines
  1. # Source Generated with Decompyle++
  2. # File: in.pyc (Python 2.4)
  3.  
  4. import math as _math
  5.  
  6. try:
  7.     (True, False)
  8. except NameError:
  9.     (True, False) = (1, 0)
  10.  
  11.  
  12. def chi2Q(x2, v, exp = _math.exp, min = min):
  13.     '''Return prob(chisq >= x2, with v degrees of freedom).
  14.  
  15.     v must be even.
  16.     '''
  17.     if not v & 1 == 0:
  18.         raise AssertionError
  19.     m = x2 / 2.0
  20.     sum = term = exp(-m)
  21.     for i in range(1, v // 2):
  22.         term *= m / i
  23.         sum += term
  24.     
  25.     return min(sum, 1.0)
  26.  
  27.  
  28. def normZ(z, sqrt2pi = _math.sqrt(2.0 * _math.pi), exp = _math.exp):
  29.     '''Return value of the unit Gaussian at z.'''
  30.     return exp(-z * z / 2.0) / sqrt2pi
  31.  
  32.  
  33. def normP(z):
  34.     '''Return area under the unit Gaussian from -inf to z.
  35.  
  36.     This is the probability that a zscore is <= z.
  37.     '''
  38.     a = abs(float(z))
  39.     if a >= 8.3000000000000007:
  40.         sum = 0.5
  41.     else:
  42.         sum2 = term = a * normZ(a)
  43.         z2 = a * a
  44.         sum = 0.0
  45.         i = 1.0
  46.         while sum != sum2:
  47.             sum = sum2
  48.             i += 2.0
  49.             term *= z2 / i
  50.             sum2 += term
  51.     if z >= 0:
  52.         result = 0.5 + sum
  53.     else:
  54.         result = 0.5 - sum
  55.     return result
  56.  
  57.  
  58. def normIQ(p, sqrt = _math.sqrt, ln = _math.log):
  59.     '''Return z such that the area under the unit Gaussian from z to +inf is p.
  60.  
  61.     Must have 0.0 <= p <= 1.0.
  62.     '''
  63.     if p <= p:
  64.         pass
  65.     elif not p <= 1.0:
  66.         raise AssertionError
  67.     flipped = False
  68.     if p > 0.5:
  69.         flipped = True
  70.         p = 1.0 - p
  71.     
  72.     if p == 0.0:
  73.         z = 8.3000000000000007
  74.     else:
  75.         t = sqrt(-2.0 * ln(p))
  76.         z = t - (2.3075299999999999 + 0.27061000000000002 * t) / (1.0 + 0.99229000000000001 * t + 0.044810000000000003 * t ** 2)
  77.     if flipped:
  78.         z = -z
  79.     
  80.     return z
  81.  
  82.  
  83. def normIP(p):
  84.     '''Return z such that the area under the unit Gaussian from -inf to z is p.
  85.  
  86.     Must have 0.0 <= p <= 1.0.
  87.     '''
  88.     z = normIQ(1.0 - p)
  89.     return z + (p - normP(z)) / normZ(z)
  90.  
  91.  
  92. def main():
  93.     Hist = Hist
  94.     import spambayes.Histogram
  95.     import sys
  96.     
  97.     class WrappedRandom:
  98.         
  99.         def __init__(self, baserandom = random.random, tabsize = 513):
  100.             self.baserandom = baserandom
  101.             self.n = tabsize
  102.             self.tab = [ baserandom() for i in range(tabsize) ]
  103.             self.next = baserandom()
  104.  
  105.         
  106.         def random(self):
  107.             result = self.next
  108.             i = int(result * self.n)
  109.             self.next = self.tab[i]
  110.             self.tab[i] = self.baserandom()
  111.             return result
  112.  
  113.  
  114.     random = WrappedRandom().random
  115.     
  116.     def judge(ps, ln = _math.log, ln2 = _math.log(2), frexp = _math.frexp):
  117.         H = S = 1.0
  118.         Hexp = Sexp = 0
  119.         for p in ps:
  120.             S *= 1.0 - p
  121.             H *= p
  122.             if S < 9.9999999999999998e-201:
  123.                 (S, e) = frexp(S)
  124.                 Sexp += e
  125.             
  126.             if H < 9.9999999999999998e-201:
  127.                 (H, e) = frexp(H)
  128.                 Hexp += e
  129.                 continue
  130.         
  131.         S = ln(S) + Sexp * ln2
  132.         H = ln(H) + Hexp * ln2
  133.         n = len(ps)
  134.         S = 1.0 - chi2Q(-2.0 * S, 2 * n)
  135.         H = 1.0 - chi2Q(-2.0 * H, 2 * n)
  136.         return (S, H, ((S - H) + 1.0) / 2.0)
  137.  
  138.     warp = 0
  139.     bias = 0.98999999999999999
  140.     if len(sys.argv) > 1:
  141.         warp = int(sys.argv[1])
  142.     
  143.     if len(sys.argv) > 2:
  144.         bias = float(sys.argv[2])
  145.     
  146.     h = Hist(20, lo = 0.0, hi = 1.0)
  147.     s = Hist(20, lo = 0.0, hi = 1.0)
  148.     score = Hist(20, lo = 0.0, hi = 1.0)
  149.     for i in range(5000):
  150.         ps = [ random() for j in range(50) ]
  151.         (s1, h1, score1) = judge(ps + [
  152.             bias] * warp)
  153.         s.add(s1)
  154.         h.add(h1)
  155.         score.add(score1)
  156.     
  157.     print 'Result for random vectors of 50 probs, +', warp, 'forced to', bias
  158.     print 
  159.     print 'H',
  160.     h.display()
  161.     print 
  162.     print 'S',
  163.     s.display()
  164.     print 
  165.     print '(S-H+1)/2',
  166.     score.display()
  167.  
  168.  
  169. def showscore(ps, ln = _math.log, ln2 = _math.log(2), frexp = _math.frexp):
  170.     H = S = 1.0
  171.     Hexp = Sexp = 0
  172.     for p in ps:
  173.         S *= 1.0 - p
  174.         H *= p
  175.         if S < 9.9999999999999998e-201:
  176.             (S, e) = frexp(S)
  177.             Sexp += e
  178.         
  179.         if H < 9.9999999999999998e-201:
  180.             (H, e) = frexp(H)
  181.             Hexp += e
  182.             continue
  183.     
  184.     S = ln(S) + Sexp * ln2
  185.     H = ln(H) + Hexp * ln2
  186.     n = len(ps)
  187.     probS = chi2Q(-2 * S, 2 * n)
  188.     probH = chi2Q(-2 * H, 2 * n)
  189.     print 'P(chisq >= %10g | v=%3d) = %10g' % (-2 * S, 2 * n, probS)
  190.     print 'P(chisq >= %10g | v=%3d) = %10g' % (-2 * H, 2 * n, probH)
  191.     S = 1.0 - probS
  192.     H = 1.0 - probH
  193.     score = ((S - H) + 1.0) / 2.0
  194.     print 'spam prob', S
  195.     print ' ham prob', H
  196.     print '(S-H+1)/2', score
  197.  
  198. if __name__ == '__main__':
  199.     import random
  200.     main()
  201.  
  202.